%function [xw BIC]=prewhitenPH(x,p,q); 
%
%Fits an ARMA model and reports estimated residuals
%
%INPUT:
%
%x time series
%p number of autoregressive coefficients 
%q number of moving-average coefficients (default 0)
%
%If no p and q are provided, the time-difference is reported
%
%
%OUTPUT:
%xw whitened time series
%BIC Bayesian information criteria is optionally computed and returned

%P. Huybers
%Harvard, 2021

function [xw coeff MdlCov]=prewhitenPH(x,p,q); 

if nargin<3,
    q=0;
end;
if nargin<2,
    p=0;
    xw=diff(x); 
    return;
end;

if p+q==0,
    xw=x;
    coeff=NaN;
    return;
end;

Mdl2=arima(p,0,q);
try, 
    [MdlEst MdlCov]=estimate(Mdl2,x(:),'Display','off');
    xw=infer(MdlEst,x(:));
    if nargout>1,
        s=summarize(MdlEst); 
        coeff.BIC=s.BIC;
        dummy=s.Table.Value;
        coeff.p=dummy(2:p+1); %autoregressive coefficients
        coeff.q=dummy(p+2:p+q+1); %moving average coefficients
        coeff.v=dummy(end); %variance
    end;
catch,
    xw=zeros(size(x))*NaN;
    coeff.BIC=NaN;
    coeff.p=ones(1,p)*NaN;
    coeff.q=ones(1,q)*NaN;
    coeff.v=NaN;    
end;

